Vaccination and variants: Retrospective model for the evolution of Covid-19 in Italy

The last year of Covid-19 pandemic has been characterized by the continuous chase between the vaccination campaign and the appearance of new variants that puts further obstacles to the possibility of eradicating the virus and returning to normality in a short period. In the present paper we develop a deterministic compartmental model to describe the evolution of the Covid-19 in Italy as a combined effect of vaccination campaign, new variant spreading and mobility restrictions. Particular attention is given to the mechanism of waning immunity, appropriately timed with respect to the effective progress of the vaccination campaign in Italy. We perform a retrospective analysis in order to explore the role that different mechanisms, such as behavioral changes, variation of the population mobility, seasonal variability of the virus infectivity, and spreading of new variants have had in shaping the epidemiological curve. We find that, in the large time window considered, the most relevant mechanism is the seasonal variation in the stability of the virus, followed by the awareness mechanism, that induces individuals to increase/relax self-protective measures when the number of active cases increases/decreases. The appearance of the Delta variant and the mobility variations have had instead only marginal effects. In absence of vaccines the emerging scenario would have been dramatic with a percentage difference in the number of total infections and total deaths, in both cases, larger than fifty per cent. The model also predicts the appearance of a more contagious variant (the Omicron variant) and its becoming dominant in January 2022.


Introduction
The Covid-19 pandemic dramatically impacted on all aspects of world population life, causing a global lasting damage at economic, social and educational level and an enormous loss in terms of human lives. The extraordinary effort made worldwide for the development of vaccines against Covid-19 allowed an equally extraordinary achievement, such as the approval of the first vaccines in one year since the beginning of the pandemic. Such a circumstance made concrete the possibility of finally defeating the virus, limiting further use of emergency measures, such as lock-downs, no longer economically sustainable. The first vaccine authorized by the US Food and Drug Administration for distribution in the United States and by the European Medicine Agency (EMA) for the European Union (EU) countries were the mRNA-vaccine Comirnaty (BNT162b2), produced by BioNTech/Pfizer (December 2020) and soon after the one produced by Moderna (mRNA-1273). At the same time, the adenovirus viral vector vaccine Russian Sputnik V, produced by the Gamaleya Research Institute of Epidemiology and Microbiology, was authorized and distributed in different countries. One month later another adenovirus viral vector vaccine, the Vaxzevria (ChAdOx1-S), developed by the Oxford University and produced by Astrazeneca was authorized by EMA for distribution in EU. Afterwards, the single-dose Janssen COVID-19 vaccine was allowed to be distributed in the US and soon after (in March 2021) its distribution authorized also by EMA for the EU countries.
However, the vaccination campaign suffered a series of setbacks, followed by successive accelerations, due to a number of circumstances. Firstly, the suspensions-and consequent restrictions of use [1,2]-of the Astrazeneca and Janssen vaccines, due to rare cases of unusual blood clots with low blood platelets occurred in some vaccinated subjects [3], slowed down the vaccination campaign during the administration of the first doses in Italy. The summer holiday period was characterized by a significant reduction of the number of doses/per day administrated. It should be added that a communication not always clear and effective by the competent bodies, with the succession of different and sometimes antithetical recommendations of use for specific population age groups, also generated skepticism in a hard core of population that still shows reticence to undergo the vaccine. On the other hand, the gradual bureaucratic strengthening of the green pass that took place over the last few months has certainly led many previously reticent individuals to get vaccinated.
In Italy the vaccination campaign started in January 2021, firstly with sanitary personnel, and subsequently by age groups. At the time writing (31 January 2022) the percentage of Italian population fully vaccinated (one dose for Janssen vaccine and two doses for the other ones) is 77.4%, the percentage of those that received at least one dose is 84.3% and those with the booster dose 56.4%.
Nowadays there are two relevant aspects that significantly affect the success of the vaccination campaign against Covid-19 pandemic. The first one is the effectiveness of different vaccines on emerging variants of the virus [17,18] and the second one is the mechanism of waning immunity, i.e. the decline of the vaccine efficacy as time passes [19,20]. Both overlapping mechanisms lead to the potential occurrence of breakthrough infections among vaccinated individuals.
Most of the vaccines actually in use were developed against the virus wild-type and tested on large scale on Alpha variant (lineage B.1.1.7). However, the Delta variant (lineage B.1.617.2), first detected in India in late 2020, has spread worldwide, becoming soon the dominant strain in United Kingdom (approximately 90-95% cases from 7 to 21 June 2021) [21]. The European Centre for Disease Prevention and Control [22] confirmed the Delta variant to be dominant in the EU at the end of August (99, 6% prevalence with CI 72-100%). In Italy, according to the Italian National Institute of Health (ISS), the Delta variant increased in less than 5 months from 1% to 97.7% at the end of August.
Analysis by Public Health England [21] and by EPIcx lab in France [23] estimated the Delta variant to be at least 60% more transmissible than the Alpha variant and the vaccines to be less effective (after a single dose it was observed a 14% absolute reduction in vaccine effectiveness against symptomatic disease with Delta compared to Alpha, and a smaller 10% reduction in effectiveness after two doses [24]), whereas similar vaccine effectiveness against hospitalization was seen with the two variants. Also in Israel, Delta variant became dominant in July (with 90% prevalence) [25], with data on Pfizer vaccine effectiveness against hospitalization essentially in agreement with British data, but with a 30% absolute reduction in effectiveness against disease [26]. More optimistic values for the effectiveness of vaccination against symptomatic disease caused by the Delta variants are reported in [27], where it is evaluated at 88.0% for mRNA vaccine and to 67.0% for viral vector vaccine.
Both for the Comirnaty [28] and Moderna vaccine [29,30], a booster dose with the same vaccine with the original SARS-CoV-2 spike protein has been strongly recommended to enforce the protection also against the Delta variant.
As stressed in [31] it is hard to differentiate the effectiveness reduction against new variant from the natural decay of immunity as time passes. According to the retrospective cohort study conducted in USA in the previously cited paper, the reduction in vaccine effectiveness against Covid-19 infections over time is probably primarily due to waning immunity with time rather than the Delta variant escaping vaccine protection.
In [32] the waning effectiveness of the vaccine in England, at 20 weeks or more after vaccination, was estimated to be 44.3% with Vaxzevria and 66.3% with Comirnaty, while the vaccine efficacy against hospitalization and death was confirmed. Similar results were found in data collected from the Israeli national database [33], and in Qatar [34]. It has been also ascertained that the immunity response decreases faster in elderly people than in young individuals [35].
This paper is devoted to the retrospective analysis of the evolution of the Covid-19 in Italy, keeping into account the vaccination rate, the variant spreading and the immunity decay. Our choice to focus on a retrospective analysis, rather than a previsional model, is motivated by the fact that the SARS-Cov2 virus has highly variable characteristics, with frequent mutations (almost 1 new variant appearing every 4 months), high sensibility to seasonal variation and rapid decay of the immunity, acquired both by vaccination and contagion. Furthermore the interplay of the virus characteristics with highly variable restrictive measures, individual behavioral changes, and oscillating progress in the vaccination campaign, generates strong stochastic effects, which can hardly be predicted within a mean field compartmental model. Our aim is instead to catch the concrete evolution that the epidemic would have had if some of the conditioning factors had not been active, while keeping all the others at work. Indeed the comprehension of the dominating mechanisms in the past epidemic spreading, within really occurred scenarios, may be extremely useful to inform sanitary and restrictive policies.
Our objective is thus to build a mean field model, as much realistic as possible, which is able to reproduce the experimental data with excellent agreement in the largest possible time window and, afterwards, to turn off one mechanism at a time in order to weigh the role that the specific mechanism have had in forging the evolution of the epidemiological curve.
To this purpose we generalize the model developed in [36], with the introduction of the appropriate compartments for vaccinated individuals and suitably modify the parameters in order to simulate the increase in prevalence of the Delta variant, starting from mid-May and becoming dominant by August. Other variants, different from Alpha and Delta are not considered in the present model. The effect of waning immunity as time passes since the completion of the vaccine cycle, obtained by constructing an efficacy decay function and by performing the appropriate averaging procedure among individuals, and the occurrence of breakthrough infections are also included in the model.
We will show that the most relevant mechanism is the seasonal variation in the stability of the virus, followed by the awareness mechanism, that induces individuals to increase/relax self-protective measures when the number of active cases increases/decreases. The appearance of the Delta variant and the mobility variations have had instead only marginal effects. The model also predicts the appearance of a more contagious variant (i.e. the Omicron variant) and its becoming dominant in January 2022.

Materials and methods
We use a compartment model similar to the one developed in [36] with the addition of appropriate vaccinated compartments. Furthermore, in order to follow the differences in the epidemic evolution among vaccinated (v) and unvaccinated (u) individuals, we split most of the compartments of the previous SEI A I S I D RD model in the vaccinated and unvaccinated sectors, indicated by the index i, with i = v, u. In particular individuals are divided in fifteen mutually exclusive classes according to their epidemiological status: the susceptible compartment S(t), the exposed vaccinated and unvaccinated compartments E i (t) (i.e. individuals that have been infected but are not yet infective), the vaccinated compartments V j (with j = 1, 2, 3 indicating the first dose, the full vaccination, i.e. the second dose except for the single dose Janseen vaccine, and the third dose, respectively), the asymptomatic infected compartments I i A ðtÞ, the symptomatic infective compartments I i S ðtÞ, the diagnosed compartments I i D ðtÞ, the dead compartments D i (t) and finally the recovered compartment R(t), that includes both unvaccinated and previously vaccinated individuals. We assume that healed individuals can loose immunity over time and return to the susceptible compartment, after an appropriate time. We do not consider demographic birth and (not related to the virus) death process. The model is calibrated on the official data of the Italian outbreak from February 20 to December 16, 2021, reported daily by ISS and publicly available in: https://github.com/pcmdpc/COVID-19 and https://github.com/owid/covid-19-data/tree/master/public/data/ vaccinations. The resulting ordinary differential equation (ODE) Eq (1) have been solved using the SciPy libraries with initial conditions reported in Table 1 of Appendix. The best fit parameters,

The compartmental SEVI A I S I D RDS model
The epidemic dynamic is governed by the fluxes of individuals among the compartments shown in Fig 1, and fully described by a system of fifteen coupled first order differential The system is closed and positive, i.e. all the state variables take non negative values for t � 0, if initialized at time 0 with non negative values, and satisfy the mass conservation law hence the sum of the states (the total population) is constant. In our model individuals move to the vaccinated compartments V j only when the vaccine protection becomes effective, two weeks after the inoculation.

The model parameters
Let us briefly review the main characteristic of the parameters. During the first wave of the pandemic, the world had to face with a completely novel virus and it took some time to understand its mechanism of action, effective therapies and treatment modalities. As a consequence, most of the parameters that governed the evolution of the epidemic in the first period significantly changed in time. After more than one year, the knowledge and the experience in preventing, diagnosing and treating the infection made some of these parameters to reach an almost stable value. However some epidemiological parameters, as for instance individual mobility, seasonal stability of the virus, risk perception, immunity etc., are intrinsically time dependent. In order to fix them, we follow the same approach as in [36], i.e. we try to avoid step-functions and look for reasonable functional behaviors to describe their evolution.
In the following we describe the principal time dependent parameters, while the time dependence of the remaining is discussed in Appendix and the constant ones are reported in Tables 2 and 5.
The model has two different transmission rate parameters, α and γ, which govern the transmission of the virus respectively from undiagnosed and diagnosed individuals. The transmission from undiagnosed symptomatic or asymptomatic individuals is still the dominant mechanism in the spread of the epidemic. Following [37][38][39], we assume the transmission rates from undiagnosed vaccinated individuals, α v , to be different (and smaller) with respect to the corresponding rates for unvaccinated people, α u . Widely proven isolation protocols allow to assume the transmission of the virus by diagnosed cases, γ, to be residual and equally effective both from infected vaccinated and unvaccinated individuals, thus we do not differentiate the parameter γ between u and v individuals.
The parameters α i can be both factorized in a pure contact term, α c (t), describing the probability per unit time that a susceptible individual meets an infected individual, and the susceptibility term, σ i (t), which takes into account the probability that a potentially contagious contact between a susceptible and an infected individual leads to a new infection. We assume α c to be independent on the vaccination status because, during the period under investigation, selective mobility restrictions for unvaccinated individuals were not yet at work. The α i can thus be written as Both terms in the previous equation are in principle subject to changes. Restrictive measures, such as lock-down and limitations of access to places and services, certainly affect the contact rate, α c , as it happened during the first period of the pandemic. In [36] the global effect of such  (7) of Appendix. The susceptibility terms σ i (t) depend on different factors, partly related to the behaviors of individuals-e.g. more careful use of self-protective measures such as face masks, hand washing, due to increased risk perception during the rising phases of the epidemic-and partly related to the characteristics of the virus such as seasonal variation in the stability of the virus in airborne [41][42][43] and increased transmissibility due to the appearance of new Covid-19 variants. The functional form of the susceptibility function is thus chosen as where σ aw is the term encoding the awareness mechanism due to risk perception, σ term encodes the modification of transmissibility due to seasonal variability, σ var encodes the emergency of new variants and the last factor, s i 0 , differentiates between vaccinated and unvaccinated individuals. Fixing s u 0 as a fitting parameter, we assume the transmissibility from vaccinated cases to be lower (s v 0 ¼ 0:5 � s u 0 ), according to studies on household transmission [37,44]. Recent literature seems to question the viral load of vaccinated infected individuals to be lower, however its decline seems to be faster than for unvaccinated individuals [45].
Following [36] we assume a prevalence based mechanism of rising awareness that increases the risk perception, inducing individuals to adopt more protective behaviors, with the effect of reducing the transmissibility of the virus during the peak. Thus we assume the awareness mechanism to act on the transmission rate by reducing it with a factor inversely proportional to the number of infective detected individuals, without any temporary effect of amplification or falsification [46]: The reduction of contagiousness during the warm season, due to a potential decline of the stability of the virus in warm environment [42][43][44][45], already considered in [36], seems to be further confirmed by more recent literature [47][48][49]. Thus we mimic the susceptibility decrease/ increase, in spring and autumn respectively, through appropriate hyperbolic tangent functions reported in Eq (8) of Appendix.
Finally the term σ var is the one encoding the susceptibility increase due to the circulation of new and more infective variants. In particular, following the ISS data [50] concerning the prevalence of Delta variant in Italy, we assume an increase, Eq (9), regulated in such a way to mimic the exponential transition from the susceptibility of the Alpha variant to the enhanced susceptibility of the Delta variant, that became dominant at the end of July (91.4% on July, 26) Due to the circulation of the Delta variant, also the parameter γ increases according to the evolution of σ var as in Eq (10).
Vaccinated people receive partial immunity against the infection, with an effectiveness that increases with the number of doses and decreases with the time occurred since the last inoculation and eventually with the appearing of new variants. In particular, the first dose gives only negligible immunity (estimated around 30% for both BNT162b2 and ChAdOx1 nCoV-19 vaccines in England [27]), whereas the third dose gives optimal immunization also against the Delta variant [51] (estimated around 90% for BNT162b2 mRNA Vaccine in Israel). For what concerns the second dose, as previously discussed, different timing leads to extremely different levels of protection among the double dose vaccinated individuals.
Following [31] we assume the reduction in vaccine effectiveness against Covid-19 infections to be primarily due to the waning immunity mechanism, rather than the Delta variant escaping vaccine protection. Thus we do not consider a further reduction factor due to new variant. It has been nowadays established [52] that the effectiveness of vaccines against the infection decreases over a period of 6 months from the inoculation of the second dose (with an estimated decay of 40%), whereas the efficacy against severe manifestation is subject to a minor decay (less than 15%). This makes it difficult to quantify the average vaccine protection level to be included in an average field model such as the one considered in the present work. To this purpose we introduce a susceptibility reduction function for vaccinated people in Eq (1) defined as f j = 1 − E j , where E j (with j = 1, 2, 3) is the vaccine efficacy after one dose, full vaccination and booster dose, respectively. Following literature, efficacy of first dose and booster dose are fixed equal to E 1 = 30% and E 3 = 90%, respectively. The evaluation of the mean efficacy over the fully vaccinated individuals, E 2 , is instead more complicated due to different timing of full vaccination. The waning in time of the Pfizer-BioNTech vaccine (which is the most used in Italy in the first phase of the vaccination campaign) was evaluated for the US veterans in Ref. [52]. We find that these data are fitted with good approximation by the following hyperbolic tangent function where e 1 , e 2 , t E and τ E are the parameters (shown in Tables 3 and 4), and t fv is the timing of full vaccination (in [52], the month when full vaccination is administered; here instead t fv is equal to 14 days after full vaccination). The time-depending mean efficacy of full vaccination is thus evaluated averaging over the fully vaccinated individuals (from which the number of individuals, who received the booster dose, was suitably subtracted) as follows: where N fv (t, t fv ) are individuals that received full vaccination at (t fv − 14) days and have not received yet booster dose at (t − 14) days.

Results
The ODE Eq (1) have been solved using the SciPy libraries with initial conditions reported in Table 1 of Appendix. In the following, the state variables, normalized over the whole population, are initialized to the values observed on Italian population on February 19, 2021. The same quantities, i.e. the fraction over the whole population, are evaluated for Italian data as reported in Table 1. The best fit parameters, obtained minimizing the χ-square with respect to the experimental data, are listed in Tables 2-4 of Appendix. Besides the detected cases, our model predicts a relevant number of undetected symptomatic and asymptomatic cases, as shown in Fig 2 (lower right). According to our prevision, the percentage of undetected cases evolves from a minimum value of 3.6% at the beginning of summer, when the outbreak slowed down, to a value larger than 43% at the end of autumn. This circumstance can be related to the tendency within the large pool of vaccinated individuals to avoid testing when asymptomatic.
This mechanism is also useful to understand the result on the relative incidence among unvaccinated and vaccinated individuals. As shown in Fig 3, our model predicts a different impact of the epidemic on vaccinated and unvaccinated individuals, in good agreement with data published by the ISS. In the inset of Fig 3, the simulated relative incidence among unvaccinated detected cases and vaccinated detected ones (blue line in figure) is compared with the same quantity evaluated on the total cases (both detected and undetected) (green line), and with the experimental incidence (red dots), clearly evaluated only on detected cases. The blue line is typically slightly higher than the red dots, and systematically higher than the green line. This last circumstance is consistent with the hypothesis that, differently from unvaccinated individuals that are frequently required to test, vaccinated individuals, specially if asymptomatic, may remain undetected more frequently than unvaccinated ones.
When comparing model predictions and epidemiological data, it must be taken into account that both are not error-free. The error on real epidemiological data is difficult to assess, due to the stochasticity inherent in the epidemic spreading (an effect completely neglected in a deterministic compartmental model in a finite population with short-range

PLOS ONE
Vaccination and variants: Retrospective model for the evolution of Covid-19 in Italy interactions), but also due to the method of data acquisition that can introduce random, but also systematic, errors (see for instance the discontinuity present in the experimental data in mid-June due to an incorrect communication of the healed individuals by some Italian regions). A source of error is certainly linked to the diagnostic capacity, which has changed significantly during the epidemic and which is mainly linked to the amount of swabs that can be done daily.
Concerning the simulations, the main source of error is the estimation of epidemic control parameters, but also of the initial conditions. In order to reproduce the experimental curves, we had to fix many different parameters, concerning the efficacy of vaccine, the waning immunity mechanism, the spreading of new variants, and other parameters such as those appearing in Eqs (9) and (5). It is reasonable to wonder how much the results obtained depend on the initial conditions and on these parameters set a priori. To explore these aspects, we carried out a sensitivity analysis focusing on one crucial number, such as the total number of cases diagnosed on the last day of the simulation, and evaluating how this number changes under modifications of the external parameters and initial conditions. We find that the results obtained are very stable for a fairly wide variation of the external parameters. In particular, Fig 4 shows that the total number of detected cases at the final day of simulation (December 16, 2021) varies less than 1% by changing the vaccine efficacy parameters in Eq (5) over a reasonable range of variability. Fig 4 represents the total detected cases as a function of the maximum efficacy of full vaccination, e 1 + e 2 , and the medium time of antibody decay, t E . It shows that the same total incidence is obtained when a decline of efficacy is counterbalanced by a suitable growth of the waning time. A similar variation less than 1% is observed by fixing e 1 + e 2 and changing the decay interval, τ E , in Eq (5), together with t E (S1 Fig). Similarly it is interesting to explore what happens by varying the parameters regulating the insurgence of the Delta variant and the corresponding augmented transmissibility in Eq (9). In particular by varying the increase in transmissibility, Δ σ , in the range [0.4, 1], the time of appearance of the new variant, t var , in the range [65,105] and the time regulating the speed in the spreading of the new variant, τ var , in the range [10,50] (Fig 5), the total number of detected cases at the final day of simulation varies less than 1%. Furthermore, it is observed that for values of Δ σ small enough, it is completely irrelevant to vary t var and τ var in the established ranges. Similar findings are obtained varying the percentages of asymptomatic infected individuals, � v and � u , or the initial conditions, in a reasonable range.
In conclusion, the errors on the simulated data due to the estimation of initial conditions and of the parameters regulating efficacy of vaccine, waning immunity, spreading of the Delta variant and fraction of asymptomatic individuals, in realistic ranges of variation, do not seem to be greater than a few percent.
As previously seen, there are many different mechanisms that contribute to the evolution of the epidemic, such as the mobility of the population, the perceived risk, the seasonal variation in the virus stability, the appearance of new variants and, obviously, the vaccination campaign.
All the previous mechanisms are encoded in the model, through appropriate terms. It is thus interesting to analyze the contribution of each term to the epidemic evolution. Fig 6 shows the trend of diagnosed cases in different scenarios obtained turning off one term at a time in Eqs (2) and (3) (and only for the σ var also in Eq (10)), or assuming the absence of vaccine. With the word "scenario" we do not intend a previsional hypothetical evolution, but a concrete evolution that the epidemic would have had if some of the conditioning factors had not been active, while keeping all the others at work. In particular, in Scenario I we disregard the effect of seasonal reduction of virus stability and infectivity during the summer period, in Scenario II we 'freeze' the awareness of individuals to the initial value, when the risk For each scenario we evaluate the increase in the total infected detected cases and in total deaths with respect to the experimental initial values of ISS at time t 0 (DI scen ¼ I scen tot ðt f Þ À I ISS tot ðt 0 Þ and DD scen ¼ D scen ðt f Þ À D ISS ðt 0 Þ) and in Table 6 we report the percentage variations of these quantities in the specific scenario with respect to those of the reference simulation (ΔI scen − ΔI ref )/ΔI ref and analogously for the death term). As already observed in [36], we find that the seasonal modulation of virus transmissibility is an essential ingredient in order to explain the evolution of the epidemic curve during the summer period. Indeed by assuming this term to be constant (violet line in Fig 6), and equal to its value in winter time, an agreement between data and simulation would have been obtained only in the first peak, then the curve for diagnosed cases, after the achievement of a minimum at the beginning of June, would start to rise again, reaching a second and significantly higher peak in summer, followed by a decrease to a long and still high plateau. In such a circumstance at the end of the simulation period (December 16, 2021) the number of total infections and total deaths would have been significantly higher (175.1% and 95.7% respectively) than the observed ones. Analogously, if the awareness mechanism hadn't been at work, a larger and higher peak would have been reached (blue line in Fig 6), followed by a rapid decrease to zero of diagnosed individuals already around the beginning of September, so the subsequent increase in both the contact and seasonal terms would not have been sufficient to produce an increase in diagnosed cases. Let us note that this behavior corresponds to the situation in which the epidemic had spread always with the same high level of awareness that population had at the end of February, when Italy was in lock-down. But, in correspondence of the peak, the awareness is lower than in the reference simulation (due to the form of Eq (4) and to the higher value of I D in the peak), while it is larger in the summer period. As a consequence, this Scenario, is characterized by a negative percentage variation in ΔI tot with respect to the reference simulation, while the

PLOS ONE
percentage variation in ΔD is positive due to the enlargement and rise of the first peak, in a phase of epidemic with higher mortality. Both contributions of the contact term, Eq (7), and of the Delta variant, Eq (9), have minor and very similar impacts on the epidemic evolution and their shutdown determines only the reduction of the summer peak and a slower autumn growth (brown and green lines in Fig 6,  respectively).
A separate analysis deserves the Scenario V with no vaccines. What would have happened in absence of vaccines? Fig 6 shows the time evolution of active cases in absence of vaccinations (orange curve). During the first peak the curve does not differ much from the experimental one, consistently with the fact that the number of vaccinated individuals was still contained in this time window. However, from the summer period, and specially in autumn, when the vaccination campaign reaches a significant percentage of the population, scenario V foresees a significant increase compared to the case of the reference simulation/experimental data. In particular according to our simulation, if vaccines were not available, at the end of simulation, one would have obtained a percentage increase of 63% in the total detected infections and of 55% in the total deaths, in qualitatively agreement with [53].
So far we have presented the model previsions up to mid-December. If we let the simulation run beyond this moment, the agreement between model predictions and real data becomes gradually worse. Actually this is not surprising because of the presence of the Omicron variant that became more and more important at the end of the year. On the other hand, the disagreement between the model and reality can give us a quantitative estimate of the presence of Omicron variant in Italy. In Fig 7 the discrepancy between the diagnosed active cases provided by

PLOS ONE
Vaccination and variants: Retrospective model for the evolution of Covid-19 in Italy the ISS and the simulated data, I D , are plotted as a function of time. If we attribute the difference between observed and expected data to the presence of the Omicron variant, we can conclude that it appeared in Italy around the beginning of December, grew rapidly and became dominant in mid-January. This picture is in substantial agreement with the data provided by ISS [54], which gives prevalence estimates at the national level on January 3, 2022 for the Delta variant 19.22% (range: 0.0%�66.7%) and for the Omicron variant 80.75% (range: 33.3%� 100%), despite the prevalence percentages are measured in a non-reproducible way within our model (the virus is sequenced by sample on positive swabs). The data in Fig 7 may be indeed influenced by a different permanence of the individuals affected by the two variants of the virus in the compartment of the diagnosed individuals, for example if individuals with the Delta variant had on average a more severe disease, they would remain longer than the individuals affected by Omicron variant in the I D compartment, systematically affecting the prevalence of the two variants estimated in this way.

Discussion
As discussed in previous sections, the focus of the present model is not to predict the epidemic evolution, but to achieve a better understanding of the virus dynamics through a retrospective analysis aimed at weighing the role of different contributions (such as mobility of the population, perceived risk, seasonal variation in the virus stability, appearance of new variants and vaccination campaign) in forging the epidemic curves. We indeed believe that the intrinsic stochasticity of the epidemic spreading of the SARS-Cov2 virus, due to its features (frequent mutations, high sensibility to seasonal variation, rapid decay of the immunity), makes the virus dynamics highly variable and previsions scarcely reliable in respect to other infective diseases. Thus we prefer to focus on the comprehension of the dominating mechanisms in the occurred epidemic spreading, in order to inform sanitary and restrictive policies.
There are many articles in recent literature regarding the COVID-19 epidemic evolution in Italy based on compartmental models and calibrated on the experimental data during the first wave of the epidemic [55][56][57][58][59][60][61] and later waves [53,62,63]. In the following we will discuss the differences with some papers that appear more similar to our work for methods and purposes.
In [63], authors develop an interesting multi-variant model with crossing immunity, aimed at studying the coexistence and transition among different variants of the virus. They also introduce a vaccination compartment, assuming that vaccine gives full immunity against each variant of the virus. The Italian model is simplified in order to consider two variants (wild and Alpha variants), the second one assumed to be more contagious than the first one, and becoming prevalent in April 2021. The agreement between the simulated and experimental data is very satisfactory for the 2020 waves of epidemic, but it shows discrepancies in respect to the experimental data of late winter-spring 2021, mainly attributed by the authors to scarce data availability and difficulties intrinsic to SIR type models to estimate the evolution of an epidemic for a long period.
In [53] authors propose an extension of the SIDARTHE model [61], with the introduction of the appropriate vaccination compartment and different transmission rates, corresponding to new variants. The paper proposes several scenarios for different vaccination schedules, in order to produce previsions of the epidemic evolution in Italy from April 2021 to January 2022. The model assumes that vaccinated people, which have successfully obtained immunity, persist in the vaccinated immune compartment with full immunity, without waning immunity mechanism.
Finally, in [64] authors propose a mean filed model with the inclusion of vaccination compartments, aimed at studying the effect of different vaccination strategies in France. The model proposes a simplified approach, in which the vaccine efficacy is reduced in a certain age-class (65 years and older) or in the whole population in a discontinuous way at a specific time (i.e. when the Delta variant becomes dominant in France). Following the authors, this reduction can be interpreted as a waning immunity, with the assumption that the individuals, within the specific age group or, less realistically, in the whole population, were vaccinated simultaneously, or in the second scenario as a reduced efficacy of vaccine due to the Delta variant becoming dominant.
Despite the apparent similarities with these articles, we believe that the novelty of the present work is in the attempt to treat as realistically as possible with all the elements discussed above, in particular with regard to the transition among different variants and the waning immunity mechanism. Indeed in our approach i) new variant appears in a continuous way and has the effect of increasing the susceptibility of the entire population, both vaccinated and unvaccinated; ii) all the fluxes among vaccinated compartments are deducted from the effective progress of the vaccination campaign in Italy (i.e. they are input data); iii) since waning immunity as time passes occurs regardless from the circulation of new variants, it is clearly distinguished from the effect due to the insurgence of new variants; iv) finally, since the vaccination campaign runs along the entire simulation time window (10 months), the waning immunity is implemented in a continuous way, which takes into account different timing of vaccinations among individuals, i.e. the effective time evolution of the fluxes among the vaccination compartments. In this way we are able to reproduce the experimental data with an excellent agreement in a very large time windows, to predict the appearance of a more contagious variant and to weigh the role that the different mechanisms have had in forging the evolution of the epidemiological curve.
The model has some limitations. The first one is that the system is assumed to be closed and protected from the injection of new cases from abroad. This circumstance is not fully justified, specially in the summer period, when the touristic flows increase. However, according to data published by ISTAT (National Institute of Statistics), even if the international tourist flow in 2021 was in recovery with respect to the year 2020 (+ 22, 3%), it was still far from the levels of 2019 (−38, 4%) [65]. The extension of the model to open system is left to future work.
Secondly, it has been shown that the vaccine efficacy to protect against severe infections is higher than the efficacy against mild or asymptomatic infection. Our model doesn't distinguish the symptomatic cases according to the severity of symptoms, being pauci, mild and severe symptomatic cases, as well as hospitalized cases, all included in the symptomatic compartment. It would be interesting to rearrange the model in order to measure differences in hospitalization and severity of symptoms between vaccinated and unvaccinated individuals. However this would involve the introduction of new compartments into the model, circumstance that we avoided in order to keep the model simpler in the present paper.
In [66] authors study the attenuation of antibody titres after the second dose, showing that the most important factor in determining the waning immunity is sex, age and smoking. Our model does not take into account the age structure of population and not even the sex groups, and thus it is not able to capture differences in infections, mortality and recovering among different groups of individuals. We leave this interesting in-depth analysis to a future work.
Finally, it will be interesting to extend the stochastic models proposed in Refs [67,68] for the pandemic H1N1 to the forth wave of Covid-19 epidemic: in this case indeed the almost total absence of mobility restrictions, which were instead present during the previous waves, would allow to apply the social contact hypothesis [69] in order to reproduce the epidemic spreading.

Conclusion
The objectives of the present work is to achieve a detailed comprehension of the main mechanisms that have been dominant in the epidemic spreading of SARS-Cov2, to be used to inform sanitary and restrictive policies and the vaccination campaign. The research questions are: What would have happened if there were no mobility restrictions? What if the virus wouldn't have exhibited such a seasonal variability? What if there were no vaccine available? Which are the most conditioning factors in shaping the epidemic spreading? In the present paper we proposed an epidemic mean field model with the introduction of vaccinated compartments, with different number of doses, and many different conditioning factors (virus mutations, sanitary polices and restrictive measures, behavioral changes of individuals and vaccinations campaign) in order to better understand the weight of different mechanisms and to answer to the above questions. One of the peculiarities of our model is to include a realistic mechanism of waning immunity that, to our knowledge, is completely new in literature.
The model obtained in this way is able to reproduce the epidemic spreading in Italy during the third and in part the fourth wave of Covid-19 with excellent agreements. The analysis of the ingredients that must be taken into account in order to reproduce the epidemiological curves teaches a lot about the disease. The strong seasonal trend of the epidemic has been confirmed, together with the role of awareness mechanisms that allow to mitigate the epidemic spreading through individual protective behaviors adopted when the risk perception increases. The effects of the appearance of the Delta variant and the contact increase during the summer months were instead only marginal, causing a slight rise in the summer peak, being the mitigating effect of summer temperatures stronger. The model also predicts with remarkable accuracy the appearance of the Omicron variant and its becoming dominant in January 2022. According to our model, in absence of a vaccination campaign, the total number of infections and deaths would have been dramatically higher, confirming the fundamental role of vaccines in containing the pandemic and saving human lives.

Appendix-Further parameters and initial conditions
In the following we describe the parameter time dependencies, not already discussed in Section 2.2, while the constant ones are reported in Tables 2 and 5. It should be stressed that most of the epidemiological compartmental models in literature include many compartments and a relevant number of parameters that regulate fluxes between these compartments. However few available observables (i.e. data sets) in general do not allow to univocally fix the parameter values in the phase space and the set of parameters exhibited is only one of the possible that, within the specific model/formulation, allow to reproduce the epidemic evolution.
Notice that for the time dependence of the parameters we preferred the use of hyperbolic tangents (except for Eqs (9) and (11), reported below, in which we adopted the exponential law, obtained by fitting the Italian data for the prevalence of Delta variant). However, while in Eq (5), the hyperbolic tangent comes up as a fit of the vaccine efficacy data, published in Ref. [52], in all other cases, the hyperbolic tangent was preferred because it allows to describe the crossover between two different values of the parameter by means of a continuous function, with derivatives of any order that are continuous in turn.
The hyperbolic tangent transition function for the pure contact term, α c , is given by where the function is modulated in such a way to produce a doubling of the average number of contacts among individuals between February/March 2021 and the autumn period (i.e. c 1 + c 2 = 2 and c 1 − c 2 = 1), and t c and τ c are fitting parameters (the values of parameters are shown in Tables 3 and 4). The susceptibility decrease/increase, in spring and autumn respectively is described by: where t term 1 , t term 2 , t term 1 and t term 2 are fitting parameters and t � is any time such that t term 1 ⪡ t � ⋘ t term 2 , such that the two hyperbolic tangent functions assume the same asymptotic value. For simplicity we put t � ¼ ðt term 1 þ t term 2 Þ=2. The variant circulation factor in the susceptibility function is modelled as Following [21, 23], Δ σ is fixed a priori equal to 0.60. The γ parameter increase due to Delta variant is expressed as: where γ 0 is a fitting parameter. Further parameters are listed below.
• χ j (with j = 1, 2, 3), the vaccination rates corresponding to first, second (or single dose for the Janseen vaccine), and booster doses. As discussed in the introduction these parameters changed discontinuously during the last year, thus we fix them as step functions through the best fits of the experimental data. In Fig 8 the simulated evolution and the time series of vaccinated individuals, as reported by ISS, are compared.
• δ, the inverse mean latent period assumed to be the same both for vaccinated and unvaccinated people. However, the emergency of new variants has been typically characterized by a reduction of the latent and incubation period [70]. Therefore we modulate its value according to the dominant variant: where Δ 1 and Δ 2 are fitting parameters, and t var and τ var are the same parameters present in Eq (9).
• y i A , y i S and y i D , with i = v, u, the recovery rates respectively for asymptomatic, symptomatic, diagnosed vaccinated and unvaccinated individuals. The recovery rates for asymptomatic and symptomatic undiagnosed individuals are assumed not dependent on time. Those of diagnosed individuals have a significant variability over time, being affected by different factors, such as the test rate (increasing when the daily number of tests increases) and the number of active cases (decreasing when the sanitary system is overload). We assume y u D and y v D to be functions of time, not depending on the vaccination status of individual, i.e.
y u D ¼ y v D ¼ y D , and following the form: y D ¼ y 1 þ y 2 tanh ððt À t y Þ=t y Þ; ð12Þ where θ 1 , θ 2 , t θ and τ θ are fitting parameters. This behavior corresponds to a decreasing in time of the number of days spent by infected diagnosed individuals in the diagnosed compartment, due to increased testing efficiency of the healthcare system through the involvement of analysis centers and pharmacies in testing operations.
• k i D and κ i , with i = v, u, the mortality rates for diagnosed and undiagnosed infected individuals, respectively. The mortality rate for unvaccinated diagnosed individuals is assumed to follow the form: where κ 1 , κ 2 , t mor 1 , t mor 2 , t mor 1 and t mor 2 are fitting parameters andt is any time such that t mor 1 ⪡t ⪡ t mor 2 , such that the two hyperbolic tangent functions assume the same asymptotic value. For simplicity we can putt ¼ ðt mor 1 þ t mor 2 Þ=2. This behavior corresponds to a

PLOS ONE
Vaccination and variants: Retrospective model for the evolution of Covid-19 in Italy decreasing of the mortality during the summer period (from κ 1 + κ 2 to κ 1 − κ 2 ), probably due to the lower age of infected people, which returns to the same value of winter/spring in autumn.
Following the data of ISS, the mortality rate for vaccinated diagnosed individuals is assumed to be lower than for unvaccinated people. In particular we assume the k v D ¼ k u D =11:5, according to the relative risk estimation among vaccinated and unvaccinated individuals reported in [71]. According to [71], the relative mortality risk between vaccinated and unvaccinated individuals is 9.0, 11.5 and 30.3, respectively, for individuals fully vaccinated more than 120 days before, for those fully vaccinated less then 120 days before and for those with booster dose. During the period under investigation, the number of individuals with booster dose was negligible and mostly of the population recently completed the vaccination cycle. Finally the mortality rate, κ i , for both vaccinated and unvaccinated undiagnosed individuals, is assumed equal to zero; • � i , with i = v, u, corresponding respectively to the asymptomatic percentages of infections among vaccinated and unvaccinated individuals. Although there is evidence that vaccinated individuals are protected from severe illness, to our knowledge, a systematic study of differences in the asymptomatic fraction of infection between vaccinated and unvaccinated individuals is still lacking. For this reason, we fix the fraction of asymptomatic vaccinated and unvaccinated individuals to be equal.
• Z i A and Z i S , with i = v, u, corresponding respectively to the detection rates for asymptomatic and symptomatic infective individuals. These parameters are chosen to be constant in time and dependent on the vaccination status of individuals, both larger for unvaccinated than vaccinated people, reflecting the tendency of unvaccinated individuals to swab more easily than vaccinated ones.
• ϕ, corresponding to the rate at which recovered individuals loose immunity. On the time scales here considered, the presence of healed individuals, who return to being susceptible with a rate between 180 −1 and 270 −1 day −1 (as estimated in literature), turns out to be totally irrelevant, so we choose to set it equal to zero.
Fitting rates present in Eq (1), whose value is not dependent on time, are given in Table 2. Constants present in Eqs (7)-(13), obtained as fitting parameters, are listed in Table 3. Days and intervals of time present in Eqs (7)-(13), obtained as fitting parameters are reported in Table 4 (in the simulations the time 0 corresponds to the initial time, t 0 reported in Table 1, i.e. February 19, 2021). Further parameters are reported in Table 5. Fig 8 shows the comparison between the simulated evolution and the time series of vaccinated individuals, as reported by ISS.
The values of R 0 , D 0 , V 10 , V 20 , V 30 , reported in Table 1, correspond to the current values at the beginning of the period under examination, t 0 ; since the number of vaccinated individuals is negligible at t 0 , I v D 0, I v A 0, I v S 0 and E v 0 are put equal to zero; the initial values for I u A 0, I u S 0 and E u 0 , which are not experimentally observables, are obtained extrapolating simulation results from previous model developed in [36]; finally S 0 is obtained as difference between 1 and the sum of the initial values of all the other compartments.